Improvement of the estimation of the infiltration function in surface irrigation systems

Surface irrigation systems are widely used on the North China Plain. The design of surface irrigation systems can be improved by developing simulation models including the advanced trajectory, recession trajectory, and infiltration time. Therefore, the objectives of this study were as follows: (1) to evaluate different models to simulate the advanced and recession trajectories, (2) to propose a new method that reduces the required observation data for estimating the infiltration time, and (3) to evaluate the accuracy of the proposed infiltration function based on the modified infiltration time function. Field experiments were conducted. The results indicated that the power function can represent the advanced and recession trajectories well. A modified function that describes the infiltration time has a high correlation and accuracy with the measured data and can be used to estimate the infiltration time. The proposed infiltration function based on the modified infiltration time function is accurate and can be used to estimate the infiltration function.


Introduction
Evaluating surface irrigation systems mainly depends on the accuracy of estimating the infiltration function [1][2][3][4][5][6].The infiltration function can be estimated based on the advance time, recession time, and infiltration time.Water advance time is the relationship between time and distance in waterfront advance, and the water recession time is the relationship between time and distance that begins when the depth of surface water upstream decreases to zero and continues until the surface is drained [7][8][9][10][11].The infiltration time is the duration it takes for the water to reach a certain depth in the soil.
Several researchers have used the power function to estimate the advanced curve [4,[6][7][8][12][13][14][15][16] and the recession curve [7,8].The parameters of the power function can be estimated using different methods: in the two-point method, the parameters can be adjusted from observations at two points of the advance/recession curve, [4,13] best matches the advanced and recession data performed by least squares [12,16], the regression for the observed advanced and recession data [7,8], and as a part of a full mathematical model using the numerical solution where the value of the parameters was not stated, but they were estimated at each measuring station as part of the numerical solution [15].Nie, Li [17] used a secondorder equation to estimate the recession time and obtained the equation's parameters through multipoint regression fitting using observed recession data.
In some studies, the value of the empirical exponent of the power function for the advance curve was determined using furrows of different slopes [18], different types of soils [19], and various other experimental conditions [20][21][22][23].In other studies, the power function's parameters for advance and recession curves were calibrated based on limited data sets (i.e., one furrow or one border) without identifying the range of these parameters, and the established parameters were used to estimate the infiltration function in the remaining treatments (furrows/borders).However, the variation in these parameters between treatments should be considered and the effect of this variation on the estimated infiltration function should be identified.Furthermore, other functions can also represent the advanced time and recession time, which may be more accurate.
The infiltration time is a critical variable used to estimate the infiltration function.The infiltration time was calculated using simple or complicated mathematical models [24][25][26][27][28] or simply by subtracting the advanced from the recession times [7-9, 17, 29, 30].In both methods, the accuracy of the infiltration time depends on the accuracy of the estimated advance and recession times and therefore on the measurement of the advance and recession times.Measuring advanced times is easier than measuring recession times, especially in a large-scale project where collecting data requires more field researchers [9].As a result, the estimated advanced time is more accurate than the estimated recession time.Therefore, reducing the required data and testing other functions to estimate the advanced time and recession time could increase the accuracy of the estimated infiltration function.Hence, the objectives of this study were as follows: (1) to evaluate different mathematical models to estimate the advanced time, recession time, and infiltration time, (2) to improve a mathematical model that increases the accuracy and reduces the required observation data to estimate the infiltration time, and (3) to evaluate the accuracy of the proposed infiltration function based on the modified infiltration time function.
This paper is organized as follows.Section 2 presents different functions to estimate advance, recession, and infiltration times, the suggested method to optimize the empirical parameters of those functions, and the suggested method to estimate the infiltration function.Then, the experimental data used to illustrate and validate the previous methods are described.Section 3 includes the results of the function's parameters of the advance, recession, and infiltration times and a new method for estimating the infiltration function is proposed.In Section 4, the results of the function's parameters are provided and the application of the proposed method to estimate the infiltration function is demonstrated.Section 5 summarizes the main conclusions.

Advanced time
The water advance curve is formulated as an empirical power equation [8,14,17,19,[31][32][33] as follows: where t adv is the time of advance (min); x is the distance of advance (m); and a and b are the empirical coefficients, which can be determined for each experiment by minimizing the differences between the observed and estimated (Eq (1)) advance times as follows: where S adv(i) and O adv(i) are the simulated (Eq (1)) and observed advance time at time i, and N is the number of data pairs.

Recession time
The recession time is approximated by an empirical power equation [7,8], an exponential equation, or a second-order parabolic equation [31] as follows: (1) Power equation (2) Exponent equation (3) A second-order parabolic equation where t rec is the time of recession (min); x is the distance of advance (m); T is the total time of the advance, recession, and depletion (min) at the head of fields; and c, d, m, n, s, j, and w are the empirical coefficients, which can be determined for each experiment by minimizing the differences between observed and estimated (Eq (3), Eq (4), or Eq (5)) recession times as follows: where S rec(i) and O rec(i) are the simulated (Eq (3), Eq (4), or Eq (5)) and observed recession time at time i, and N is the number of data pairs.

Infiltration times
The infiltration times (Δt, min) can be obtained based on the differences between advanced and recession times.Hence, the infiltration times can be determined as a function of x by subtracting Eq (1) from Eq (3), Eq (4), or Eq (5).In the Results section, the recession equations (Eq (3), Eq (4), and Eq (5)) are compared, and the equation that gives the best results is used to obtain the infiltration times.
The infiltration times can also be approximated by a second-order parabolic equation as follows: where Δt is the infiltration time (min); x is the distance of advance (m); and f, g, and h are the empirical coefficients, which can be determined for each experiment by minimizing the differences between the observed and estimated values (Eq (7)) infiltration times as follows: where S Δt(i) and O Δt(i) are the simulated (Eq (7)) and observed infiltration time at time i, and N is the number of data pairs.

Infiltration function
The infiltration function is described by the Kostiakov equation, which can be written as follows: where Z is the cumulative infiltration depth, Δt is the opportunity time k, and α are empirical coefficients (k> 0, and 0 � α<1).
The simple postirrigation volume balance (PIVB) [34,35] is used to solve the k and α parameters based on the following equations [36][37][38]: Another way to solve the Kostiakov equation using PIVB is by defining the k and α parameters as a function of t 100 , the time needed to infiltrate 100 mm [36].

Border irrigation experiments
In this study, border irrigation experiments were conducted in the town of Nanpi in Hebei Province, China, at a longitude, latitude, and elevation of 116˚40'E, 38˚06'N and 20 m, respectively.Winter wheat was planted in October 2014 and harvested in the middle of June 2015.
The soil texture at the site is classified as silt loam (67.02% silt, 25.19% sand, and 7.79% clay on average), with a mean bulk density of 1.39 g/cm 3 for the 1 m soil depth.

Infiltration parameters
The Kostiakov equation parameters obtained by the WinSRFR model [10,11,39] based on Eq (10), and Eq (11) can simulate the border irrigation process well.The input data included the irrigation requirement, the advanced and recession times at each measurement distance within each border (every 5 m along the border length), the border dimensions (measured using a measuring tape in the field), and the elevations along the border length (measured using an optical level) were recorded for border irrigation.The average reported slope in each block varied from border to border and throughout the irrigation season.Salahou, Jiao [40] published the details of the calculation process for a closed-end border irrigation performance evaluation model considering a distance-based cutoff, and the calculation results are only summarized in this section.Since the infiltration parameters obtained were well calibrated, the distribution of the infiltration water amount simulated by the WinSRFR model along the border length can be regarded as the measured results, which can be used to test the estimated infiltration water amount distribution function in this study.The details of the layout of the border and measuring points can be found in [40].The details of the border irrigation are listed in Table 1

Statistical analysis of the results
Two statistical indicators were used to analyze the goodness-of-fit between the observations and simulations: (1) the root mean square error (RMSE) and ( 2) the coefficient of determination, which are calculated as follows: RMSE ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi 1 N where S i represents the predicted value, O i represents the observed value, � O is the mean value of all observed events at time i, and N is the number of data pairs.The RMSE indicates the agreement between the predicted and observed values.R 2 ranged between 0 and 1.
In addition to the previous statistical indicators, maximum, minimum, mean, standard deviation (σ), and coefficient of variation (Cv) [41] were used to describe the estimated parameters.X m , σ, and Cv can be calculated using the following equations: where σ is the standard deviation; xi is the estimated parameters of the ith point; X m is the mean value; N is the number of data points; and C v is the coefficient of variation, which reflects the variation in the estimated parameters (i.e., when C V < 0.1 represents weak variability; 0.1 � CV � 1.0 represents moderate variability, and C V > 1.0 represents strong variability) [42].A statistical analysis was conducted on the parameters of the power function (a, and b) in Table 2.The maximum, minimum, mean, range, standard deviation σ, and coefficient of variation Cv are shown in Table 2.The variation range of coefficient a was 0.0096-0.0535min�m -b , the mean was 0.0278 min�m -b , and the coefficient of variation was 0.370, which indicated moderate variability.The variation range of coefficient b was 1.4465 to 1.7186, the mean is 1.5495, and the coefficient of variation was 0.043, indicating relatively weak variability.

Advancement time
Fig 2 shows the advance curve of borders B02, B06, B11, and B21.The fitting curve exhibited a high degree of agreement.Even after the water cut off, the water advancement process was well in line with the power function relationship.The results showed that the advance process conformed to the power function relationship.The R 2 was between 0.9961 and 0.9997, indicating a very strong correlation.The RMSE did not exceed 0.98 min, and the curve fitting was excellent.

Recession time
The power function (Eq (3)) was used to fit the surface water flow recession process in border irrigation.The exponential function (Eq (4)) was also used to fit the recession process of the border irrigation events.The variation range of the empirical parameters m and n is small, and both increase with the decrease in inflow rate (Fig 4).
The third method for estimating the recession curve involves a second-order parabolic function (Eq (5)).Fig 5 shows the variation range of the empirical parameters s, j, and w.The parameter s was small and equal to zero at most borders, converting Eq (5) into a straight-line equation.[13] represented the recession curve by a straight-line equation, but they divided the trajectory into two straight lines.
The statistical characteristics of the parameters of the power function (Eq (3)), the exponential function (Eq (4)), and the second-order parabolic function (Eq (5)) are shown in Table 3.Among the parameters, the variation range of the border recession time T was the largest, ranging from 18.83 to 55.50 min, and the second was the coefficient m, which ranged from 25.075 to 57.680 min.The coefficients s and n have the smallest variation range, and the ranges were only 0.001 and 0.006, respectively.In terms of variability, the coefficients of variation of the coefficients c and s were the largest, at 0.569 and 1.505, respectively, indicating strong variability, and the other coefficients all showed moderate variability.especially at the head of the border.In addition, due to the limitation of the formula itself, the exponential function cannot describe the recession process of the convex shape of some borders well.The R 2 by the exponential function (Eq (4)) was between 0.147 and 0.956.Except for B07, the R 2 was low, only 0.147, and the correlation was not significant.The slope of the first part of the B07 border (0-15 m) was steep, while the first half (15-50 m) was low-lying, and the water infiltration time of the field surface here was long, which did not conform to the trend of exponential function change, resulting in a low fitting correlation coefficient.The RMSE was between 1.31 and 5.52 min, which did not differ greatly from the RMSE in the result of the power function (Eq (3)) (2.53 min).The power function fitting curve did not have a problem describing the shape of the recession curve.The R 2 by the power function was between 0.605 and 0.9997, showing a very significant correlation.The RMSE was between 0.14 and 4.28 min, which was slightly larger than the fitting result of the advanced process.The second-order parabolic function (Eq (5)) described the recession curve well.The R 2 was between 0.155 and 0.0.957,showing a very significant correlation.The RMSE was between 1.31 and 5.33 min, which was similar to that of the exponential function.
Therefore, compared with the previous three methods (Eq (3), Eq (4), and (Eq (5)), the exponential function (Eq (4)) has a simpler form and does not require the determination of the beginning of the recession phase (T).However, the formula form of the exponential function (Eq (4)) determines that the fitting curve of the recession process can only be a downward convex shape, which is inconsistent with the measured data in some fields.The results of 16 borders showed that the power function was better than the other functions.The average RMSE fitted by the power function was 2.53 min lower than that of the other functions (Eq (4), and (Eq (5)), and the average R 2 was 0.858 higher than that of the other functions.

Infiltration time
Based on the previous results, Eq (3) best represents the recession trajectory.Hence, we determine the infiltration time as a function of x by subtracting Eq (1) from Eq (3) as follows: We noted in the previous sections that b and d are not strong empirical parameters, and the difference between their averages is approximately 0.8, so we can rewrite Eq (19) as follows: The infiltration time is described by Eq (7), Eq (19), and Eq (20), respectively, and compared with the measured results.The RMSE and R 2 results are shown in Table 4.Among the above three methods, the RMSE conforms to the following relationship: method 1 > method 2 > method 3, and the correlation coefficient shows the opposite law.The calculation results of Eq (20) have the best correlation and the highest accuracy with the measured data, which can describe the infiltration time curve of border irrigation.This method (Eq (20)) is used to estimate the infiltration function in the following section.

Infiltration function
The soil infiltration process can be described by the Kostiakov model Eq (9), and the infiltration time at each point in the border field satisfies Eq (20).Hence, the distribution formula of the infiltration water amount along the border length direction can be expressed as follows: As mentioned in the methodology section, the distribution of the infiltration water amount simulated by the WinSRFR model along the border length can be regarded as the measured results, which can be used to test the estimated infiltration water amount distribution function using Eq (21).The infiltration functions in the 16 borders are described by Eq (21), and the parameters (k, and α) of Eq (21) were calculated based on Eq (13), and Eq (12) (See section 2.4  for obtaining the Kostiakov parameter, and see section 3.3 for fitting coefficient values such as T, p, q, and r) Table 5 shows the results of the RMSE and R 2 values of the infiltration function (Eq (21)) of the 16 borders.The mean square error RMSE is between 3.32 and 6.94 min, and the average value of R 2 is 0.825.

Formula
The infiltration functions of the B02, B06, B11, and B21 borders are shown in Fig 7 .The calculated value of the infiltration function using Eq (21) is in good agreement with the simulated value using the WinSRFR model with little difference.Hence, we conclude that the infiltration water distribution function established in this work is reasonable, and the calculation results are accurate.

Discussion
The advance power law with an exponent (b) between 1.44 and 1.72 can represent a different advance rate with time, with an average of 1.55, which is within the basic range of the previous studies.Amer and Amer [7], Amer [8] estimated the advanced power law with an exponent of 1.42.The range of b values reported in other studies under various conditions.For example, Serralheiro [18] reported b values between 0.52 and 0.94 for Mediterranean soil for furrows of different slopes.Alvarez [19] estimated b values between 0.58 and 0.72 for different types of soil for furrow irrigation systems.In another study, Khatri and Smith [20] reported b values for 27 furrows with an average of 0.86, and Ebrahimian and Liaghat [22] reported b values for 5 furrows and 6 borders with averages of 0.63 and 0.58, respectively.Therefore, the b value is a function of several factors, such as soil moisture, soil type, surface slope, and irrigation system.Hence, fitting power functions to field data increases the accuracy of the estimated infiltration function.
Relationships for estimating recession times in blocked-end borders have been previously proposed [35,43], but they can only be used with prior knowledge of infiltration.This process consists of two distinct phases.In the upper part of the field, water will be lost primarily through free drainage, and recession will occur relatively quickly.For the lower portion of the field, water can only be lost by infiltration.Hence, the recession will take much longer and will be approximately linear with time if the slope is uniform.In the real world, recession times will depart from linearity due to the nonuniformity of the bottom slope.comparisons of measured and simulated recession time values based on the different fitting procedures tested.The measured data are somewhat atypical of blocked-end borders, in which some ponding occurs at the end of the field.In these tests, the inflow was cut off and carefully timed to prevent ponding [40].A fundamental problem with the PIVB method is that one equation is available to solve for all unknown parameters of the infiltration equation being used.Even when dealing with the two parameters of the Kostiakov equation, an infinite number of solutions can be proposed.To overcome this problem, an empirical relationship developed by Merriam and Keller [34] was used.In this study, the developed empirical relationship for recession times is applicable only to the conditions of the study.To apply this method under other conditions, recession times must be measured, and a new empirical relationship must be developed for those new conditions.Thus, the problem of needing to measure recession times in the field was partially solved.The recession time can be measured in some experiments, and then the method presented in this study can be applied to the remaining experiments.Although the recession time measurements are inaccurate in comparison in advance, which, in turn, justifies the proposal to use estimates instead of measured values, the advance time measurements are subject to considerable uncertainty just as recession ones, even though those measurements are more precise.The arrival of the irrigation stream to a measuring station is affected by other factors besides infiltration, e.g., hydraulic resistance, bottom slope undulations, soil clods, vegetation, cracking, etc.In a border, the stream typically does not advance uniformly across the width.Thus, advanced measurements may be precise, but they may not accurately represent infiltration.While recession time measurements are imprecise, they provide more information about long-term infiltration rates than the advance times.The former only provides information about the early part of the infiltration process.
Recession times can be difficult to obtain manually, but with current developments in electronics, we can expect to develop sensors that not only make it easier and less costly to collect recession data but also measure flow depths during the postadvance phase of irrigation.Those depths are critical to improving the accuracy of infiltration estimates and to determining hydraulic resistance independently from infiltration.However, sensors are not available at many research sites.

Conclusion
Establishing the infiltration function in a surface irrigation system is critical.Different methods were presented to estimate the advanced trajectory, recession trajectory, infiltration time, and infiltration functions.The results showed that the advanced trajectory can be represented well with the power function Eq (1).By comparing the three methods for estimating the recession trajectory, the power equation Eq (3) best fit the observation data.When the total time of the advance, recession, and depletion at the head of borders was not collected, the secondorder parabolic equation Eq (5) can be an alternative to estimate the recession trajectory.The results showed that the s parameter of the second-order parabolic equation was small and equal to zero at most borders, converting Eq (5) into a straight-line equation.The modified function (Eq (20)), which describes the infiltration time of border irrigation, had the best correlation and the highest accuracy with the measured data and, therefore, can be used to estimate the infiltration time.A new method was improved to estimate the water distribution along the field.The results show that the new infiltration function (Eq (21)) is in good agreement with the simulation results obtained from the WinSRFR model.Hence, the infiltration water distribution function established in this work is reasonable, and the calculation results are accurate.The border irrigation system was tested using the established equations in this study.However, a furrow irrigation system could also be tested.

Fig 1
Fig 1 shows the advance equation parameters a and b (Eq (1)).A statistical analysis was conducted on the parameters of the power function (a, and b) in Table2.The maximum, minimum, mean, range, standard deviation σ, and coefficient of variation Cv are shown in Table2.The variation range of coefficient a was 0.0096-0.0535min�m -b , the mean was 0.0278 min�m -b , and the coefficient of variation was 0.370, which indicated moderate variability.The variation range of coefficient b was 1.4465 to 1.7186, the mean is 1.5495, and the coefficient of variation was 0.043, indicating relatively weak variability.Fig2shows the advance curve of borders B02, B06, B11, and B21.The fitting curve exhibited a high degree of agreement.Even after the water cut off, the water advancement process was well in line with the power function relationship.The results showed that the advance process conformed to the power function relationship.The R 2 was between 0.9961 and 0.9997, indicating a very strong correlation.The RMSE did not exceed 0.98 min, and the curve fitting was excellent.

Fig 3
shows the results of the power function's empirical parameters (c and d).The parameters c and d (Eq (3)) were scattered, and the trend in the change with the influence of the inflow rate was not clear.The d value was less than 1.0, and the recession curve was convex.

Fig 6
shows the fitting results of the recession curve of the B02, B06, B11, and B21 borders.The exponential function fit poorly on the gradual recession process of border irrigation,